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Abstract 



The cubic nonlinear Schrodinger equation with a lattice potential is used to model a pe- 
riodic dilute gas Bose-Einstein condensate. Both two- and three-dimensional condensates are 
• considered, for atomic species with either repulsive or attractive interactions. A family of exact 

£S) \ solutions and corresponding potential is presented in terms of elliptic functions. The dynami- 

cal stability of these exact solutions is examined using both analytical and numerical methods. 
For condensates with repulsive atomic interactions, all stable, trivial-phase solutions are off-set 
from the zero level. For condensates with attractive atomic interactions, no stable solutions are 
found, in contrast to the one-dimensional case (§) . 

g . 

Y . 1 Introduction 

The cubic nonlinear Schrodinger equation with repulsive or attractive nonlinearity and an external 



q | potential models the mean- field dynamics of a dilute-gas Bose Einstein condensate (BEC) [|i9| , 28]. 

In this context, the equation is also known as the Gross-Pitaevskii equation. BECs are of interest in 
the theoretical and experimental physics community since they are examples of macroscopic quan- 



tum phenomena which manifest phase coherence M, 11, EQ, UM. BECs have possible applications in 



such different areas as quantum logic [pj, matter- wave diffraction [25], and matter- wave transport 



After rescaling of the physical variables, the governing equation is 

i t) = ~Atj;(x, t) + a \iP{x, t)\ 2 rP(S, t) + V{x)^{x, t), (1) 

where A denotes the Laplacian operator, ip(x, t) is the macroscopic wave function of the condensate 
in one, two or three dimensions, with x defined y) or (x, y, z) respectively. The real- 

valued function V(x) is an experimentally generated macroscopic potential. The parameter a 
determines whether (p]) is repulsive (a = 1, defocusing nonlinearity), or attractive, (a = — 1, 
focusing nonlinearity). The dimensionless variables in (ffl) are related to the physical variables by 
t -> (4irh\a\N/Mn)t, x -» {Q./^\a\N) 1 / 2 x, ip -> (M^)" 1 /^, andU(f) -* (Mn/An\a\N)V(x). Here 
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Table 1: Physical parameters for various BEC experiments along with the lifetime of the condensate 
and the elapsed scaled time At for ([!]). Note that both 7 Li and 85 Rb have a negative s-wave 
scattering length which imply an attractive condensate. The last column of the table gives the 
reference for the corresponding experiment. For the data in the table: * denotes the total time 
of the experiment, f is estimated from the peak width of Fig. 3 of [p2]| , £ is calculated using the 
experimental results of [3C]. 



N is the number of atoms in the volume f2 containing all atoms (condensed and uncondensed) in 
the the experiment (typically on the order of (10/iin) 3 ||), M is the mass of a single atom, and a is 
the s-wave scattering length for collisions between atoms. Depending on the atomic species, a can 
be either positive or negative so that both signs of a = sign(a) are relevant for BEC applications. 
In Table |], various physical quantities for BEC experiments are displayed. Note that the s-wave 
scattering lengths for 7 Li and 8 5Rb are negative, so these condensates are attractive, whereas the 
other experiments listed use repulsive condensates. 

Current BEC experiments require a confining potential in order to trap the condensate. For 
theoretical convenience, this confinement is often assumed to be harmonic. Once confined, addi- 
tional electromagnetic potentials may be imposed on the BEC. In particular, there has been recent 
interest in confinement of repulsive BECs using standing light waves which results in a sinusoidal 
potential in the Nonlinear Schrodinger equation ([l]). Although experiments focus on the quasi- 
one-dimensional regime, quasi-two-dimensional and three-dimensional lattice potentials have been 



suggested [23|. The interaction of the BECs with such potentials is of great practical and theoreti- 



cal interest. The quasi-one-dimensional regime [10, 27] applies when the longitudinal dimension of 



the condensate is far greater than the transverse dimensions, which are themselves of the order of 
the healing length of the condensate. Likewise, the quasi-two-dimensional regime [^] is significant 
when two dimensions of the condensate are far greater than the third dimension, which is itself of 
the order of the healing length of the condensate. The healing length is the width of the boundary 
region over which the probability density of the condensate drops to zero. Thus these quasi-one- or 
two-dimensional approximations allow for simplification of dTJ). In reality, however, all condensates 
are three-dimensional and in many cases no reduction in dimension is justifiable. 

In one dimension, the potential generated in experiments is a standing light wave Q and 
V(x) = —Vq sin 2 (mx) = (Vq cos(2m:z;) — Vq)/2. The same potential was considered in [||, 13]. 



Typically, the condensate is distributed over tens of periods of the standing light wave |1|. In 
|^. f|, [lO|| , a number of exact solutions of (Uj) with this potential was constructed and their stability 
was investigated. Additionally, a more general potential was considered: 

V(x) = -V Q sn 2 (mx,k), (2) 

where sn(x, k) denotes the Jacobian elliptic sine function, with elliptic modulus k. It is periodic 
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Figure 1: External potential V{x) which exhibits both a periodic and a confining structure. The 
periodic structure is given by sn 2 (x, k) which reduces to a standing light wave (as k — ► 0) in the 
top figure. Other possibilities of the potential (|2|) include well-separated troughs or peaks. In this 
paper, we consider the periodic region away from the edge of the confining experimental potential. 
Typically this region extends over tens of periods Q. 



with period 4K(k) = 4f£ /2 da/y/l - k 2 sin 2 a. This family of potentials contains the standing 
light wave potential as a special case: limfc_>o sn(x, k) = sin(x). This more general potential is 
considered because it allows the inclusion of different regimes of BEC dynamics. In particular, for 
values of k up to k = 0.9, the shape of the potential is virtually indistinguishable from a sinusoidal 
one, but for values of k close to 1, i.e., k > 0.999, (|2|) gives periodic potentials with well-separated 
peaks or troughs. Figure |] illustrates the different regimes of the periodic potential @. The figure 
also indicates the presence of a confining potential which is applied in all current experiments. 

Due to technical complications, there are currently no experiments where a BEC is trapped 
in a higher-than-one-dimensional periodic potential. However, the interest in the applications 
mentioned above strongly suggests that these experiments will eventually take place. Already, 
theoretical investigations of BECs in multidimensional lattice potentials suggest the realization of 
such experiments |2^]. Motivated by the previously referenced developments in BECs, we consider 
(|l]) with repulsive and attractive nonlinearity and periodic potential in two and three dimensions. 
Note that in [|16| , the repulsive condensate in two-dimensions was investigated. 

A judicious choice of the potential allows the construction of a large class of exact solutions, as 
in one dimension p, |j| [h], [l6|]. The family of potentials considered is 

d d 

V(x) = - JJ (Aisn 2 (miXi, ki) + Bi) + m 2 k 2 sn 2 (miXi, ki), (3) 

i=l i=l 

where d is the number of dimensions: d = 2 or d = 3, and x\ = x, X2 = y, and X3 = z. Here Ai, 
Bi and mi, i = 1, ...,d are real constants. The elliptic moduli ki are in the interval [0, 1]. The 
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Figure 2: Various lattice potentials. For all figures, A\ = A 2 = B\ = B 2 = 1. For (a), k\ = k 2 = 0, 
mi = m 2 = 1; for (b), fci = 0.999, fc 2 = 0, mi = 2K(0.999)/vr, m 2 = 1. Finally, for (c) 
k 1 = k 2 = 0.999 and m 1 = m 2 = 2K(0.999)/?r. 



first term is a straightforward generalization of the one-dimensional potential (|2|), with an additive 
constant. The other terms facilitate the construction of exact solutions. 

Note that in the important physical case of a trigonometric potential fcj — ► 0, i = l,...,d, 
the potential reduces to a product of standing light wave potentials. This generalization of the 
one-dimensional case is completely analogous with the standard approach of higher-dimensional 
Kronig- Penney potentials |l8[| . For non-trigonometric potentials, the exact expression (||) for the 
potential is necessary to allow the construction of exact solutions. Nevertheless, it is the qualitative 
features of the potential, i.e., its periodicity and amplitude, that are most important for the stability 
properties of the solution. As numerical and analytical results throughout this paper demonstrate, 
the behavior of a solution in a lattice potential is largely independent of the quantitative features 
of the potential, just as in the one-dimensional case || ^, |n|. Figures || and || display the potential 
(|3|) for various values of its parameters in two and three dimensions respectively. 

Because of the present lack of experiments in two or three dimensions, this paper can either 
be read as considering (|l]) with a periodic potential as a mathematical problem in its own right, 
or as a collection of predictions from the mathematical mean-field model concerning the dynamics 
and stability of BECs in higher-dimensional periodic potentials. The paper is arranged as follows: 
exact solutions of (||) with periodic potential (||) are constructed. A linear stability analysis is 
presented for these exact solutions. However, since in many cases the analysis is inconclusive, a 
complete determination of the stability properties of all our exact solutions is obtained by numerical 
simulation. We conclude the paper with a brief summary of the results and their implications for 
BEC dynamics. 

2 Exact solutions 

For typical potentials of (|l]), the construction of exact solutions is not obvious and usually not 
possible. For the potentials (g) and (||), it is possible to find families of exact stationary solutions. 
We look for solutions whose spatial dependence is separated: 

d 

iP(x,t) = e- i " t H<t> i (xi)- (4) 

4=1 

Here d = 2 or d = 3. The derivation for d = 1 is similar, but slightly easier. Inserting the ansatz 
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Figure 3: Various 3-D lattice potentials. For all figures, A\ = A 2 = B\ = B 2 = — 1. For (a), 
k\ = k 2 = k% = and mi = ra2 = 7713 = 1; for (b), k\ = k 2 = 0, £13 = 0.999, mi = m,2 = 1, 
m 3 = 2if (0.999) /V, for (c), h = 0, k 2 = h = 0.999, mi = 1, m 2 = m 3 = 2if (0.999)/vr, for (d), 
k 1 = k 2 = k 3 = 0.999, mi = m 2 = m 3 = 2K (0.999) /vr. 



@ into (0) gives 

w = 4Ex +«n^| 2 + y(x), (5) 
i=l ^ i=l 

where prime denotes differentiation with respect to the argument. In order to separate this equation, 
the potential is written as 

d d 

V(x) = -aH\4>i\ 2 + Y,W l (x i ), (6) 
i=\ i=i 

which results in the d separated equations 

Ui<t>i = -\$ + Wi{xi)^i, (7) 
with Ya=i = w - Using an amplitude-phase decomposition 

4>i{x i ) = r i {x i )e ie ^\ (8) 



eliminating the overall exponential factor, and equating the imaginary part to zero, gives 



( x =°i I jTTZK^ ( 9 ) 



dw 
/ rf{w) 

where q is an integration constant. Multiplying the remaining equation by r-, results in 

„2 



■\i r 'if - A + 2 / Wiix^nixyiix^dxi. (10) 
2 2rf J 

This equation is rewritten using the new function Si = rf. Combining this with a substitution on 
the integral term finally gives 

{S[f = -SuiSf - 4c? + 8St J Wi{Si)dSi. (11) 
Provided the right-hand side of this equation is a polynomial in Si of degree 3 or 4, the equation 



is solvable in terms of elliptic functions [36|. This imposes conditions on the potential for this 
construction to work: Wj as a function of Si is a polynomial of degree at most 2 |36|| . Using this 
result, one finds 

rf = Aisn 2 (m,iXi,ki) + B { , (12) 
which gives rise to the potentials (0). The parameters uj, A4, Bi, c, are constrained by 



E"i = E(^ ? (i+^)+^f). (13) 



i=l i=l 



c| = ^g^ + ^^ + fcf^). (14) 

For d = 2, this results in a one-parameter family of solutions: specifying the potential fixes the 
parameters A1A2, B\/A\, and -E^/^b- Thus A\jAi is a free parameter. Similarly, for d = 3, the 
cross ratio (vli : .A 2 : A3) is free, resulting in a two-parameter family of solutions. The existence of 
these solutions requires rf > and cf > 0. This imposes additional constraints on the parameters: 
Ai >0, Bi> 0, or ^ > 0, -Ai < B { < -Ai/kf, i = l,...,d. 
Note that all our solutions satisfy 

r 2K(ki) d 

iP(x + 2K(k i )e i ,t) = 4,(x,t)e^, li = c i ir r, (15) 

Jo r i( w > 

where el is the unit vector in the X{ direction. Thus the solutions are nonlinear generalizations 
of the Bloch states of solid state theory Note, however, that any such analogy is flawed in 
principle: the reasoning leading to Bloch functions in solid state theory is based on methods for 
linear differential equations, whereas the governing equation (|l]) is nonlinear. Similarly flawed is 
any interpretation of uj as an eigenvalue or a separation constant. However, for small amplitude 
solutions for which the linear terms dominate the nonlinearity, the intuition gained from solid state 
band-gap theory can be useful in interpreting the weakly nonlinear results. 

From (12), it follows that the amplitude Ti{xi) is periodic in Xi with period 2K(ki) /rrij. However, 



generally the phase 9i(xi + 2K (k{) /rrii) 7^ 6i(xi) + 2nir, for some integer n. Thus, the solution ip(x, t) 
is usually not periodic in any of the Xi, see (|i~5|). Imposing periodicity in the ^-direction requires 
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a quantization of the phase 9i{xi). It is unclear if such a quantization in all directions is possible, 
since not enough free parameters are available to satisfy the number of quantization conditions. 
There are two special cases in which phase quantization is not a concern. The first case results 
in trivial-phase solutions, for which q = 0. The second case is the trigonometric limit, in which 
ki = 0, i = 1, . . . , d. 

The solution has trivial phase in the ^-direction if q is zero. There are three possibilities: 

Bi = 0: rj = \f~Ki sn(m,iX, ki), (16a) 

Bi = -Ai : n = yJ-Ai cn(rriiX,ki), (16b) 

Bi = -Ai/kf : n = \j-Aijkl dn{mix,ki), (16c) 

where cn(mjXj,fcj) is the Jacobian elliptic cosine function, and dn(m,iXi,ki) denotes the third 
Jacobian elliptic function. 

The solution is trigonometric in the z-direction if ki is zero. Then 



r 2 (x) = Ai sin 2 (xj) + B { , tan(^) = J 1 + ^ tan(mjXj), (17) 

and phase quantization is satisfied. Notice that it is possible for the solution to have trivial phase 
in one direction and be trigonometric in the other. In contrast, as k{ — > 1 it is possible to obtain 
solitary wave solutions. These are discussed in |l0|| . More details on the trigonometric solutions 
and the trivial phase solutions is also given there. 

3 Linear stability 

The linear stability analysis of the solutions found in the previous sections is completely analogous 
to the one-dimensional case, discussed in || 1C]. The dimensionality of the solutions only requires 



minor alterations. 

Consider an infinitesimally small perturbation of an exact solution 



ip(x, t) = (r(x) + e<p(x, t)) exp {i9(x) — iut) , (18) 



with r(x) = r[j=i r «( x «) an d 9{x) = Yli=i@i{ x i) an d r i{ x i)i &i(xi) defined in (12) and @ respec- 
tively. Substituting this in (|]) and linearizing in e gives 

i^+LUip = --A<p + ar 2 (ip* + 2<p) + V(x)ip + -{V9) 2 <p - -<pA9 - iV9 ■ Vip, (19) 
(..) ( 2 2 2 

where ip* is the complex conjugate of ip, VO-Vip = Yli=i §I~§x~-> an< ^ C^@) 2 = V9-S79. Decomposing 
(p into its real and imaginary parts p(x,t) = u(x,t) + iv(x,t), (|19|) becomes 

!(:)-( VZ)(:)- 

wit h ./ j x o ) anC ^ 

L+ = -^A + 3ar 2 + V{x) -uj + ^(W9) 2 , (21a) 
L_ = -^A + ar 2 + V(x) -uj + ^(\79) 2 , (21b) 
S = -V9-V-^A9. (21c) 
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For solutions with non-trivial phase, the analysis of this linear system is currently an open prob- 
lem. We restrict ourselves to the stability analysis of trivial-phase solutions: V# = 0. Separating 
the time dependence from the spatial variations, 



D-d 



e A i, (22) 



we obtain for trivial-phase solutions 



Hv)-'ro *-°)(v> <-> 

withL + = — ^A + 3ar 2 + V(x) — uj and L_ = — ^A + ar 2 + V(x) — uj. Note that the perturbations on 
the trivial-phase solutions considered do not necessarily have a trivial-phase profile, since lm((p) = v 
is not necessarily zero. 

Adding the potential V{x) to the Nonlinear Schrodinger equation destroys three of its four 
continuous symmetries. However, equation (|l|) is still phase invariant: the transformation ^{x, t) — > 
e* 7 ^(x, t), with 7 a real constant, leaves ([!]) unchanged. Using Noether's theorem, 

L_r(x) = 0. (24) 

Thus, A = is in the spectrum of the Schrodinger operator L_: G <r(L_). Since both L + and L_ 
are self-adjoint operators, their spectrum is contained on the real axis. Furthermore, since — is 
a positive operator, and r(x), V(x) are bounded, these spectra are contained in some semi-infinite 
interval: 

a(L + ) G [A +> oo), a(LJ) G [A_, oo), (25) 

with A± = inf |W |=i <T7|^ ± |77>, ||?7|| 2 = /„ \ V \ 2 dx, Q = {(x,y) G [0,2K(fcx)] x [0,2K(k 2 )]} for d = 2, 
and n = {(x,y,z) G [0,2if(*!i)] x [0,2K(k 2 )\ x [0, 2K (k 3 )]} for d = 3. Note that since G cr(L-), 
A_ < 0. 

Theorem 1 For condensates with repulsive atomic interactions (a = 1) an exact trivial-phase 
solution such that r(x) > is linearly stable. 

Proof In this case a = 1 and L + = L_ + 2r 2 . Thus, A + > A_. If r(x) > 0, then it is the 
ground state of L_ [14] and A_ = 0, A+ > 0. Note that for one-dimensional Schrodinger operators, 
the ground state is unique. This is not necessarily true for the higher-dimensional case we are 

1/2 

considering. Since L + is a positive operator, we can construct L! . Consider 



H = L^L-L 1 / 2 . (26) 



Let f = L+ /2 [/, then 



(F + A 2 )£ = 0. (27) 
The operator H is also self-adjoint. Let //o be its smallest eigenvalue, then 

Mo = inf (v\H\ri) 
IMI=i 



inf ( 77 

Nl=i \ 

inf (l 1 / 2 !! 
\\n\\=l \ + 



i}[ 2 l^lT 



i] 



L. 



Li, J] 
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]/2 — 1 /2 / 1/2 

which implies /io > 0. Note that by letting r/ = L + r/\\L + r\\, \L^_ rj 
fxo = 0. Thus A 2 < 0, and A is imaginary or zero. 



Li, 1] 



so that 



Remark As in [1C], we have proven the nonexistence of exponentially unstable eigenvalues. 
However, algebraic instabilities of the zero or other modes are not ruled out. 

Theorem 2 For condensates with attractive atomic interactions (a = —1) an exact trivial-phase 
solution such that r(x) > is linearly unstable. 

Proof In this case a = — 1 and L_ = L + + 2r 2 . Thus, A_ > A + . If r(x) > 0, then, as above, it is 
the ground state of L_ [14] and A_ = 0, A+ < 0. Since L_ is a positive operator, we can construct 
L 1 / 2 . Consider 

H = L l ' 2 L^ 2 . 



Let £ = L V2 C7, then 

(F + A 2 )£ = 0. 

The operator H is self-adjoint. Again, let /j,q be its smallest eigenvalue, then 



(28) 
(29) 



inf (r,\H\ 
Ml =1 



inf ( r? 

Hl=i 



inf ( L}l 2 ri 

hll=i 



L 1/2 L+L 1/2 



L. 



rl/2 

L_ rj 



< 0. 

Thus there is at least one A 2 such that A 2 > 0, resulting in a pair of opposite, real A's. 



For some other cases, more can be said about the stability or instability of an exact solution, as 
in [^0|. However, those stability criteria are not as straightforward to apply as the two given here, 
and we only make use of the two criteria above. 

It follows from our derivation of (£4j) that the solution which has a dn(miXi, fcj) function in all 
directions is a deformation of the ground state of the linear Schrodinger equation with periodic 
potential given by (||). Since (||) is nonlinear, the concept of a ground state as an eigenfunction 
corresponding to the lowest eigenvalue as used in quantum mechanics has to be abandoned. How- 
ever, the ground state can still be defined as the global minimizer (if one exists) of the Hamiltonian 
of the equation, with the constraint HV'II 2 = C, for some constant C. This Hamiltonian is 

w) = X G w ■ vr + v{ ^ 2 + ^°i^i 4 )) ds - ( 3 °) 

It is bounded from below only for condensates with repulsive interactions. Appendix A (due to 
Jared C. Bronski) demonstrates that in this case, the dn(mjXj, fcj) solution is indeed the ground state 
of the nonlinear equation. Since the Hamiltonian (|30| ) is not bounded from below for condensates 
with attractive interactions, no such ground state exists for this case. 
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Density: l¥(x,y,t)l 2 Contour: l¥(x,y,t)l 2 Fourier modes 
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Figure 4: The two-dimensional stable repulsive evolution of the dn(mjXj, fej) solution corresponding 
to ( |l6cj ) over four periods, with k\ = &2 = 0.5, mi = m2 = 1, and A\ = A2 = —1. 



4 Condensates with repulsive interaction 

In this section, the numerical solutions of (|l]) in the repulsive regime of the nonlinear Schrodinger 
equation are considered. The initial conditions are selected from the exact solutions given by (Q) and 
(|l^-|r^). Of particular interest are the trivial-phase solutions given by ( 16a -c) and the nontrivial 



phase solutions in the trigonometric limit ([17]). These solutions automatically satisfy the phase 



quantization condition for our periodic domain ||, [Kj] and correspond to a solution with a ramped 
phase profile. The numerical procedure used is a fourth-order Runge-Kutta method in time and a 
filtered pseudospectral method in space [^]. For each experiment, a small amount of white noise 
was added to the initial conditions as a perturbation to expedite the onset of any instability. 

4.1 Numerical Simulations: Square, Regular Lattice 

Our computational studies begin by considering the case of a square, regular lattice. Thus the 
spatial domains of all X{ are identical as are the elliptic moduli ki, and the generated lattice 
potential V{x) and the solution ifj(x,t) do not distinguish between the different directions X{. 
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Density: l¥(x,y,t)l 2 Contoui" l¥(x,y,t)l 2 Fourier modes 




Figure 5: The two-dimensional unstable repulsive evolution of the sn(mjXj, k{) solution corre- 
sponding to (|l6a| ) over four periods, with k\ = &2 = 0.5, m\ = mi = 1, and A\ = A2 = 1. 



4.1.1 Two-Dimensional Solutions 

The first two-dimensional solution we consider is the dn(rriiXi, ki) solution (|l6cj ). As predicted by 
Theorem [I], this nodeless solution is linearly stable for a condensate with repulsive interactions. The 
evolution of four periods of this solution with Ai = —1, m; = 1 and k = 0.5, i = 1, 2 for t £ [0, 60] 
is shown in Fig. ||. The three columns of this figure, and all other figures of this type, represent 
in order the evolution of the density \ip(x, y, t)\ 2 , a contour plot of this same evolution, and the 
evolution of the arctan of the Fourier spectrum. The bottom picture in the first two columns shows 
the potential V(x,y). The arctan is applied to the spectral evolution to limit the range of the 
power spectrum. This results in a more elucidating representation of the dynamics of the exact 
solution, by removing the extreme contrasts in the range of the Fourier spectrum. For the nodeless 
dn(rriiXi,ki) solution, we note that over the time range considered, and indeed for times t — > 00, 



the solution does not change even when strongly perturbed. As argued previously [||, 10, |16|] , the 
stability of this solution is reminiscent of the stability of the plane wave solution of the nonlinear 
Schrodinger equation with repulsive nonlinearity |3^|. 



In contrast to the stable dn(rriiXi,ki) solution, the sn(mjXj,fci) and cn(mjXj,/cj) solutions cor- 



responding to ( 16a ) and ( 16b ) are unstable. Both these solutions have nodes and violate the linear 
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Density: IT(x,y,t)l 2 Contour: IT(x,y,t)l 2 Fourier modes 




Figure 6: The two-dimensional unstable repulsive evolution of the cn(mjXj,fcj) solution corre- 
sponding to (|l6b|) over four periods, with fci = = 0.5, mi = 777,2 = 1, and A\ = A2 = —1. 



stability criterion of Theorem [T|. The evolution for t G [0, 60] of four periods of these solutions is 
shown in Figs. || and ^. The evolution of the density clearly shows the onset of a modulational 
instability which deforms the exact solution. The power spectrum shows that this modulational 
instability results in the activation of a large number of Fourier modes, destroying the possibility 
of stable evolution. 

To illustrate the computational stability results, we calculate the difference E between the den- 
sity of the exact solution and the density of the numerical (perturbed) solution for the sn(m,iXi, ki), 
cn(mjXj, h L ) and dn(m,iXi, ki) solutions of ( |16a| )-( |l6c| ). Thus, E = \ip(x,y,t)\ 2 — \ip(x, y, 0)| 2 . If 
the numerical evolution was exact and the initial condition was unperturbed by white noise, this 
difference would be identically zero, since the solutions we consider are stationary. Thus any growth 
in this difference is due to an instability mechanism which causes any errors to grow. In Fig. 0, the 
three columns given by E(Dn, Dn), E(Sn, Sn), and E(Cn,Cn) represent the errors in the exact 
solutions given by the dn(mjXi, ki), sn(mjXj, ki) and cn(mjXj, ki) solutions of (|16a| )- (|16c| ) respec- 
tively. For the dn(mjXj,/cj) solution, the error E(Dn, Dn) remains at the level of the initial noise. 
In contrast, the errors for the sn(rriiXi,ki) and cn(mjXj, ki) solutions, E(Sn,Sn) and E(Cn,Cn), 
start to grow at the onset of instability near t « 40 and t ps 20 respectively. 
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E(Sn,Sn) E(Cn,Cn) 




t=60 



t=40 
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Figure 7: The two-dimensional repulsive evolution of the error, E 



(x,y,t)\< 



for 



the sn(miXi,ki), cn(rriiXi,ki) and dn(mjXj, fcj) solutions of (|16a| )-( |16q ). The solutions considered 
correspond to Figs. 0-0. The bottom row shows the respective exact solutions. 



The unstable behavior is further illustrated by the evolution of the L°° norm of E: 
maxs xy \ ||^(x,y,i)| 2 — |^(x,y,0)| 2 |. In Fig. ^, this L°° norm dynamics is given for both the 
two-dimensional and three-dimensional solutions of the equation with repulsive nonlinearity con- 
sidered numerically in this paper. The L°° norm of the unstable solutions grows at the onset of the 
instability and saturates at a finite value. This reflects the nature of the repulsive instability, i.e., 
no large gradients or sharp peaks are allowed to develop in the solution. For the stable solutions, 
this L°° norm remains at the initial noise level. 



4.1.2 Three-Dimensional Solutions 



As in the two-dimensional case, we first consider the trivial phase dn(mjXj,A;j) solution (16c). As 
shown in Theorem |l|, this nodeless solution is linearly stable. The evolution of eight periods of 
this solution with Ai = —0.5, m, = 1 and k = 0.5, i = 1,2,3 for t € [0,40] is shown in Fig. |]. 
The two rows of this figure, and all other figures of this type, represent the evolution of a selected 
density iso-surface \ip(x, y, z, t)\ 2 = constant, and the evolution of an iso-surface of the arctan of 
the Fourier spectrum. Unless otherwise stated, the density iso-surface shown is at 30% of the value 
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L (t) sn(x.k) 

L"(t) nontrivial phase stable 

L°°(t) nontrivial phase unstable 



40 80 20 40 

Time (t) Time (t) 

Figure 8: Evolution of the L°° norm of the error in two- and three-dimensions for several trivial- 
phase solutions of the repulsive equation. 



from minimum to maximum while the iso-surface of the Fourier spectrum is at a value of one. 
Note that the arctan of the spectrum limits the value of the power density to a maximum of n/2. 
Over the time range considered, and indeed for times t — > oo, the nodeless dn(mj2j, hi) solution 
is unaffected by the perturbations. This even holds for perturbations that are strong compared to 
the amplitude of the solution. 

Just as in the two-dimensional case, the sn(mjXi, k{) and cn(mjXj, hi) solutions corresponding to 
(|l6aj-b) are unstable. Both these solutions have nodes and violate the stability criterion of Theorem 
[l]. The evolution for t € [0,40] of eight periods of these solutions is shown in Figs. 10 and 11. The 
evolution of the density clearly shows the onset of a modulational instability which causes break-up 
of the exact solution. The Fourier spectrum of the density shows that this modulational instability 
results in the activation of a large number of Fourier modes, destroying the possibility of stable 
evolution. Fig. || illustrates this unstable behavior by showing the evolution of the L°° norm of 
the error: max; Ii!/)Z j | \ip(x, y, z, t)\ 2 — \i/j(x, y, z, 0)| 2 1 . As with the two-dimensional results, the L°° 
norm of the error of the unstable solutions grows at the onset of the instability and saturates at a 
finite value. For the stable solutions, this L°° norm remains at the initial noise level. 
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Figure 9: The three-dimensional stable repulsive evolution of the dn(rriiXi, ki) solution correspond- 
ing to (16c) over eight periods (two in each direction), with k\ = ki = k% = 0.5, m\ = mi = = 1, 
and Ai = A<i = A% = —0.5. The Fourier spectrum is composed of six symmetric peaks around the 
origin which are obscured by the grid lines. 



t=0 



t-20 



L=40 






Figure 10: The three-dimensional unstable repulsive evolution of the sn(mjXj,fcj) solution cor- 
responding to ( |16a| ) over eight periods (two in each direction), with k% = ki = k% = 0.5, 
mi = 1712 = m 3 = 1] and A\ = A-i = A3 = 1. At t = 40 the further loss of structure, which 
is exemplified by the Fourier spectrum, is hidden by the shape of the iso-surfaces at the boundary 
of the box. 
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4.2 Numerical Simulations: Rectangular, Irregular Lattice 

In addition to regular lattices, we can consider more complicated solutions by modifying the elliptic 
modulus ki and periodicity parameter m*. We restrict ourselves to stable two-dimensional dynamics 
given by the dn(mjXj, ki) solution of (|16c|) since they are easy to illustrate and stable solutions are 
the most relevant for experiments. The three-dimensional behavior follows in a straightforward 
manner from a generalization of these two-dimensional visualizations. We only show the stable 
solutions at the initial time since they are unchanged as t — > oo. Figure [l^ displays the various 
behaviors as the elliptic modulus is varied from k = 0.1 to k = 0.999 for various periods and values 
of ttt. j . The different stable solutions vary from periodic lattice solutions in Fig. Ilia and |i~2|c, to the 



well-separated and localized spikes of 12b, to well-separated regions of oscillation and localization 



12d-f. Thus stability is independent of the lattice structure, as long as the solution is off-set from 



the zero level. 

4.3 Numerical Simulations: Nontrivial Phase 

Analytical results for the stability of the nontrivial phase solutions are difficult to obtain even in 
one dimension. Thus we rely on numerical investigations of the stability of these solutions. To 
avoid the complications which arise from phase quantization, consider the solutions (|i~7| ) that are 
trigonometric in all directions and for which phase quantization is automatically satisfied. For 
the NLS equation with repulsive nonlinearity (jlj) in two- and three-dimensions, these solutions 
are stable or unstable depending on the offset parameters B\ and Bi (two-dimensions) or Bi,B% 



and £?3 (three-dimensions). As shown in Fig. 13 the off-set solution in two-dimensions, which is 



qualitatively like the stable dn(mjXj, ki) solution, is stable with B\ = 1 and B2 = 0.7 whereas in 
Fig. |l4| the unstable nontrivial phase solution which is below the offset threshold is illustrated with 
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Figure 12: Stable two-dimensional solutions with irregular lattice potentials and m, values chosen 
so that the solutions are 47r-periodic. (a) Two periods in each direction, with k\ = 0.1 and &2 = 0.1; 
(b) two periods in each direction, with k\ = 0.9 and k% = 0.9; (c) four periods in x and two periods 
in y, with k\ = 0.1 and k% = 0.1 respectively; (d) two periods in each direction, with k\ = 0.1 
and &2 = 0.9; (e) four periods in x and two periods in y, with k\ = 0.1 and k% = 0.9 respectively; 
(f) two periods in x and four periods in y, with k\ = 0.1 and &2 = 0.9 respectively. The spatial 
domains are normalized to x,y G [— 2ir,2ir]. 



B\ = 1 and F>2 = 0.6. These nontrivial phase solutions illustrate the necessity of offset for stability. 
In particular, for the values considered, an offset threshold is achieved at B<± ~ 0.65. 

Similar results for repulsive nontrivial phase solutions hold in three dimensions: the offset 
parameter B{ is the crucial parameter which determines the stability of a given nontrivial phase 
solution. Figure 15 shows the stability of an off-set solution in three-dimensions with B\ = B<i = 
B3 = 1 whereas in Fig. 16 an unstable nontrivial phase solution, which is below the instability 
threshold, is illustrated with B\ = B2 = -B3 = 0.5. 

Mixed solutions can also be considered: these are solutions with nontrivial phase (|l7]) in one or 
more directions and trivial phase (16a) in the remaining directions. As before, offset determines the 
stability of these solutions. To illustrate this, we consider the two-dimensional repulsive evolution 
with a trivial phase dn(mjXj, ki) solution of (|l6c|) in the y direction and a nontrivial phase solution 
(|j~7|) in the x direction. The dn(miXj,/cj) solution of (|16c| ) was found to be stable whereas the 
nontrivial phase solution ([17]) was stable provided a sufficient offset was present. Following the 
previous paragraphs, we shown in Fig. [l?] the mixed solution with offset parameter B\ = 0.4 for 
the nontrivial phase dimension. The resulting evolution is stable under perturbation as t — > 00. 
Alternatively, when the offset parameter is decreased to B\ = 0.3, the nontrivial phase solution is 
unstable, as observed in Fig. 18. 
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Figure 13: The two-dimensional stable defocusing evolution of the nontrivial phase solution ( |17| ) 
with four periods, and k% = hi = 0.0, m\ = m<i = 1, A\ = = 1, B\ = 1, and B<i = 0.7. 



5 Condensates with attractive interaction 

In this section, the numerical solutions in the attractive regime of the nonlinear Schrodinger equa- 
tion (H) are considered. The initial conditions are only selected from the exact trivial-phase solu- 
tions given by ( |16a| -c). These trivial phase solutions, given by the sn(mjXj, ki), cn(miXi, ki) and 
dn(niiXi, ki) (i = 1,2, or 3) of ( |16a[ -c), characterize the basic solution types: peak-on-peak with 
nodes, peak-on-peak without nodes, peak-on-trough with nodes, and peak-on-trough without nodes. 

It is well known |33| that in two dimensions and higher, solutions of the attractive nonlinear 
Schrodinger equation undergo a process of collapse and blow-up. This process is also possible in 
the presence of a periodic potential. It is briefly discussed in Appendix B. In addition, plane wave 
solutions are known to be modulationally unstable |33|| . Thus we do not expect stable two- and 
three-dimensional solutions. However, it may be possible to obtain a solution for which the onset 
of instability occurs on a longer time scale than the experimental lifetime of a condensate. This 
possibility is examined numerically. 
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Figure 14: The two-dimensional unstable repulsive evolution of the nontrivial phase solution ( |i"7| ) 
with four periods, with k\ = hi = 0.0, ra\ = 7712 = 1, A\ = = 1, B\ = 1, and B2 = 0.6. 



5.1 Numerical Simulations: Square, Regular Lattice 

Our numerical considerations begin with the case of a square, regular lattice. Thus the spatial 
domains of all Xi are identical as are the elliptic moduli k%. Thus, the generated lattice potential 
V(x) and the solution t/j(x,t) do not distinguish between the different directions Xj. 
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t=0 t=40 L=80 




Figure 15: The three-dimensional stable repulsive evolution of the nontrivial phase solution ( |l7|) 
with eight periods (two in each direction), and k\ = k<i = k% = 0.0, m\ = mi = = 1, 
A\ = A2 = A3 = 1, and B\ = B2 = -B3 = 1. The Fourier spectrum is composed of six symmetric 
peaks around the origin which are obscured by the grid lines. 




Figure 16: The three-dimensional unstable repulsive evolution of the nontrivial phase solution 
(|j~7|) with eight periods (two in each direction), with k\ = k2 = k% = 0.0, m\ = iri2 = = 1, 
Ai = A 2 = A 3 = 1, and B x = B 2 = B 3 = 0.5. 
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Figure 17: The two-dimensional stable repulsive evolution of a mixed solution which has nontrivial 
phase profile dl7| ) in the x direction and the dn(mjXj, k{) profile of (16c) in the y direction. The 
figure shows two periods in each direction with k\ = 0, k<i = 0.5, m\ = ni2 = 1, Ai = 1, A2 = —0.5, 
and Bi = 0.4. 



5.1.1 Two-Dimensional Solutions 

The first two-dimensional solution considered is the dn(m.jXj,/cj) solution ( |16c| ). For the attractive 
nonlinear Schrodinger equation, the dn(miXi, ki) solution is proved to be unstable in Theorem [2|. 
This is a manifestation of the modulational instability of plane waves [33|| . Although the modula- 
tional instability is present in the dynamics, the behavior is dominated by the collapse and blow-up 
phenomena which are characteristic of the attractive regime. This behavior is observed in Fig. [l^. 
In this case, ln(l + \ip(x, y, t)\ 2 ), and its contours are plotted to better illustrate the collapse and 
blow-up. The evolution of four periods of this solution with Ai = — 1, mi = 1 and k = 0.5, i = 1,2 
for t € [0, 0.56] shows how quickly this phenomenon occurs. 

In order to effectively capture the onset of collapse, we have aided the instability in some 
simulations by considering initial conditions which have twice the amplitude of an exact solution. 
Our motivation in expediting the instability is to limit computational cost: since the exact solutions 
are stationary, sufficient perturbation is required to achieve collapse. Building up such perturbations 
from our regular white-noise perturbations requires long time scales relative to the time scale on 
which collapse occurs (on the order of 1000 times longer for our numerical experiments). In order 
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Figure 18: The two-dimensional unstable repulsive evolution of a mixed solution which has a 
nontrivial phase profile (17) in the x direction and the dn(m,iXi, fcj) profile of ( |16c| ) in the y direction. 
The figure shows two periods with k\ = 0, A>2 = 0.5, mi = m,2 = 1, Ai = 1, ^2 = —0.5, and 
Si = 0.3. 



to resolve the collapse, a small time step is required, which is unnecessary for the time leading 
up to collapse. An alternative approach to avoid this problem is to use an adaptive time-stepping 
algorithm. The main effect of our approach is to alter the time of collapse and blow-up. Other 
qualitative features of the solution are unaffected. 

Note that the double-amplitude dn(mjXj, fcj) initial condition, which is nodeless and has a larger 
L 2 norm than either the sn(rriiXi,ki) or cn(rriiXi,ki) solution, collapses into a few distinct peaks. 
Additionally, the spreading of the Fourier mode spectrum reflects the nature of the localization 
inherent to collapse. 

The exact peak-on-peak sn(mjXj, fcj) solution is also unstable for attractive condensates as shown 
in Fig. 20. As with the dn(mjXj, fcj) initial condition, the solution eventually blows up near t ~ 31. 



In this case, the total L norm is much smaller than that of the dn(mjXj, fcj) solution so that only a 



single collapse peak is observed. In contrast to Fig. 19, the initial conditions are the exact solutions 
with a small amount of white noise added. 

The cn(rriiXi,ki) solution can be either peak-on-peak or peak-on-trough, depending on the 
parameters. In the absence of the mean-field nonlinearity, the peak-on-trough cn(niiXi, kj) solution 
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Figure 19: The two-dimensional unstable attractive evolution of the double-amplitude dn(mjXj, ki) 
initial condition corresponding to (16c) over four periods, with k\ = k 2 = 0.5, m\ = m 2 = 1, and 
Ai = A 2 = -1. 



would be stable. However, the cubic nonlinearity once again gives rise to collapse and blow-up. 
The evolution of a double-amplitude initial condition is shown in Fig. 21 with A\ = A 2 = — 1. 
The collapse of this evolved initial condition occurs shortly after t ~ 0.47, with the formation of 
well-defined growing peaks at the locations of the maxima of the initial condition. Since it is the 
nonlinearity which drives the collapse, we inhibit the blow-up by choosing an exact solution with a 



lower amplitude and smaller resulting nonlinearity. With A\ = A 2 = —0.5, the dynamics of Fig. |22 



is significantly different from that of Fig. 21. In particular, the onset of collapse does not occur 
until t ~ 200. Furthermore, the L 2 norm for the 16 periods considered computationally is just 
large enough to lead to collapse. Thus for small-amplitude solutions, it may be possible to obtain 
a periodic condensate in the lattice potential over the lifetime of the experiment. Note from Table 
||, that St = 200 far exceeds the time scale of any experiment with attractive condensates. 

The unstable behavior is further illustrated by the evolution of the L°° norm of the error: 
max^ y} | \ip(x, y, t)\ 2 — \ip(x, y, 0)| 2 1. In Fig. [2^, the dynamics of the L°° norm is given for both the 
two-dimensional and three-dimensional solutions of the attractive equation considered numerically 
in this paper. The L°° norm of the unstable attractive solutions grow to infinity at the onset of 
the collapse of the peaks. This illustrates the fundamental nature of the attractive collapse and 
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Figure 20: The two-dimensional unstable attractive evolution of the sn(rriiXi,ki) solution corre- 
sponding to (|16c|) over four periods, with k\ = ki = 0.5, m\ = ni2 = 1, and A\ = Ai = 1. 



blow-up, i.e., large gradients or sharp peaks develop in the solution. 



5.1.2 Three-Dimensional Solutions 

Following the discussion of the two-dimensional solutions, we consider the trivial phase dn(mjXj, hi) 
solution fll6c|) . Also in three dimensions the plane- wave solution is modulationally unstable and the 
solution is linearly unstable. The unstable behavior however is dominated by collapse and blow-up 
which are characteristic of the attractive regime. The evolution of two periods of this solution with 
Ai = -0.25, mj = 1 and k = 0.5, i = 1,2,3 for t e [0,0.114] is shown in Fig. ||. As with the 
two-dimensional solutions, we expedite the collapse and blow-up by doubling the amplitude of some 
of our exact solutions. The three rows of this figure, and of all other figures of this type, represent 
a slice in the (x,y)-plane of the evolution of the condensate density \i/j(x, y, 0, t)\ 2 , an iso-surface 
of \ip(x, y, z, t)| 2 , and the evolution of an iso-surface of the arctan of the Fourier spectrum. For the 
attractive equation, the displayed iso-surface for the density is at 50% of the value from minimum 
to maximum while the iso-surface of the Fourier spectrum is shown at a value of one. All other 
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Figure 21: The two-dimensional unstable attractive evolution of the double-amplitude cn(mjXj, kj) 
initial condition corresponding to (|16c| ) with four periods, and k\ = k 2 = 0.5, mi = rri2 = 1, and 
A 1 = A 2 = -l. 
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Figure 22: The two-dimensional unstable attractive evolution of the exact cn(m,iXi,ki) solution 



corresponding to (16c) with 16 periods, and k\ = = 0.5, rri\ = mi = 1, and A\ = A2 = —0.5. 
Note that the collapse time is well beyond t « 200. 
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Figure 23: Evolution of the L°° norm of the error in two- and three-dimensions for several trivial- 
phase solutions of the attractive equation. 
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Figure 24: The three-dimensional unstable attractive evolution of the double-amplitude 



dn(miXi, ki) initial condition corresponding to (16c) with eight periods (two in each direction), 
with k\ = &2 = &3 = 0.5, mi = ni2 = = 1, and A\ = A2 = A3 = —0.25. 
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Figure 25: The three-dimensional unstable attractive evolution of a double-amplitude sn(mjXj, ki 
initial condition corresponding to ( |16c| ) with eight periods (two in each direction), with k\ = = 
&3 = 0.5, mi = "i2 = Tn-s = 1, and A\ = A2 = A3 = 1. 
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figures of this type are similar. This figure clearly shows the localization of the collapsing solution 
which results in the spreading of the Fourier modes in the spectrum. 

A double-amplitude sn(mjXj, fcj) peak-on-peak initial condition is also unstable in the attractive 
limit as shown in Fig. [2^. As with the dn(mjXj, ki) initial condition, the dynamics eventually lead 
to blow-up near t ~ 0.07. In this case, the total I? norm is much larger than for its two-dimensional 
counterpart so that more distinct peaks are observed during the collapse process. 

As in two dimensions, an exact peak-on-trough cn(mjXj,/cj) solution is stable in the absence 
of the mean-field nonlinearity. However, the cubic nonlinearity once again gives rise to collapse 
and blow-up. The evolution of a double-amplitude initial condition is shown in Fig. ^ with 
A\ = A2 = —1. The dynamics leads to collapse of this initial condition at t ~ 0.07, and the 
formation of a few well-defined peaks is observed. Since it is the nonlinearity which drives the 
collapse phenomena, we inhibit the blow-up by choosing an exact initial condition with lower 
amplitude and smaller resulting nonlinearity. Choosing A\ = A2 = —0.5, the resulting dynamics 
of Fig. |27] is shown to be significantly different from those considered previously. In particular, 
the onset of collapse is arrested by the fact that only two periods in each direction are considered 
so that the 1? norm is small, and the solution remains in the linear regime for t £ [0,2000]. This 
again suggests the possibility of obtaining a stable periodic condensate in a lattice potential over 
the lifetime of an experiment. 



6 Discussion 

We considered the repulsive and attractive nonlinear Schrodinger equation with an elliptic function 
potential in two and three dimensions. Periodic solutions of this equation were constructed and their 
dynamical stability was investigated analytically and numerically. For condensates with repulsive 
atomic interactions, all stable, trivial-phase solutions are deformations of the ground state of the 
linear Schrodinger equation. These solutions are off-set from the zero level. This is confirmed with 
extensive numerical simulations on the governing nonlinear equation. Likewise, nontrivial-phase 
solutions are stable if they are sufficiently off-set. Thus for condensates with repulsive interactions, 
a large number of condensed atoms is sufficient to form a stable, periodic condensate, as in the 



one-dimensional case [10 [ . Physically, this implies stability of states near the Thomas-Fermi limit. 

For condensates with attractive atomic interactions, no stable solutions are found, in contrast 
to the one-dimensional case ||. However, the time scale for the onset of instability for some of our 
solutions far exceeds the time scale of current physical experiments with attractive condensates. 
This occurs when the solutions are localized in the troughs of the potential and have nodes. These 
solutions may be observable in experiments, given the time scale over which the instability develops. 

Many issues remain open. In particular, for the majority of our solutions, we have no analytical 
results concerning stability or instability. Nontrivial phase solutions especially have eluded analysis. 
A more detailed understanding of the collapse and blow up in periodic potentials is necessary. Of 



particular interest is the full class of initial conditions which lead to collapse and blow-up [31 



35]. Furthermore, the effect of a confining potential in addition to a periodic potential warrants 
investigation. These are only a sampling of the many mathematical and physical issues which arise 
in this interesting problem. 
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Figure 26: The three-dimensional unstable attractive evolution of a double-amplitude cn(mjXj, ki 
initial condition corresponding to (16c) with eight periods (two in each direction), with k% = h,2 - 
k^ = 0.5, mi = m<2 = = 1, and A\ = A2 = A3 = — 1. 
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Appendix A: The ground state of the Nonlinear Schrodinger equa- 
tion with repulsive atomic interaction 

This appendix is due to Jared C. Bronski. 

Theorem 3 A stationary solution ip(x,t) = rffle^®^ - ^ with r(x) > of (Q), a = 1 is a global 
minimizer of the Hamiltonian $3(\ ) with the constraint that \\ip(x,t)\\ 2 = C, for some constant C. 



Proof Consider the reduced Hamiltonian 

7i(V0 = W) - AHVII 2 = J • W>* + V(xM 2 + i|V| 4 - A|Vf) dx, (31) 

where A is a Lagrange multiplier. Its critical points are determined by 

SHty) = => Xtj; = --Ai/; + MV + V(x)ip. (32) 

Thus the Lagrange multiplier A, which enforces the constraint of fixed particle number, is interpreted 
physically as the chemical potential u. In what follows, the factor e~ lujt in the solution is omitted. 
In general, an expression for A = u> in terms of the critical point ip is obtained by multiplying this 
expression by ip* and integrating over Q: 

u ~ m • (33) 

Global minimizers of T~C(ip) subject to the constraint ||^>|| 2 = C may be assumed to be real, since 



\\r(x)e i6 ^\\ 2 = ||r(f)|| 2 

H(r(x)e ie ^) = H{r{x)) + ^ I r 2 (Vd) 2 dx > H(r(x)), 



(34a) 
(34b) 



with equality only valid for X79 = 0. In this case, because of the phase invariance of (JT]), we can 
choose 9 = 0. Clearly a global minimizer is only fixed up to multiplication by a constant phase 
factor. 

The second variation of the reduced Hamiltonian is given by the quadratic form 



5 2 H{ip) 



-V0 • V<f + V{x)\<j>\ z + ^V 2 + ^<t? + 2|Vf - w 



n2 ,2 



(u\L + \u) + {v\L_\v) , 



dx 



(35) 



where ip is assumed to be real, cp = u + iv, and L + and L_ are defined by ( |21a[ -b) with 8 = 0. From 
the proof of Theorem [l], it follows that for ip{x) = r{x) > 0, the operator L+ is positive definite, 
whereas L_ is positive semi-definite. Thus ip{x) = r{x) > is a local minimizer of the Hamiltonian, 
subject to the constraint \\ip\\ 2 = C. 
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The proof that a positive solution is necessarily a global minimizer is by contradiction. Assume 
there exists a ip such that \\ip\\ 2 = C and 7Y(V>) < Ti.(ip). By assumption, tp is real. This can also 
be assumed for ip since this can only lower the energy. Now consider 



^>{v) = cos{v)ip + i sin(z/)^. 
Then ^(v) satisfies the constraint ^(z^H 2 = C, since 



cos 2 (V 



+ sin 2 (z^ 



Also, 



= cos 2 (^) y Qv^-W + V(£ 



cos 2 (^) 
sin 2 (V) 
cos 2 (^) sin 2 (f 



-V$ • V^* + V(x 



C. 



2 + - cos 2 (z^ 

2 + - sin 2 (z^ 
2 v 



< cos 2 (^)W(V') +sm 2 {v)H(?p). 
This last inequality follows from f n ( cos 4 (z^)|^| 4 + 2cos 2 (V) sin 2 (z/) 1 "' 1 J 



(36) 
(37) 

dx + 
4 ^cfx + 

(38) 

j 2 + sin 4 (z/)|V;| 4 W < 



In (cos 2 (i/) 1^1 4 + sin 2 (i/)|^| 4 j dx, which is easily verified. Equation ( pg| ) expresses the convexity of 
the Hamiltonian as a functional of ip. Now, since TC(ip) < H(ip), near = 



(39) 



and thus Ti.(^f(u)) is quadratically decreasing near v = 0. But, 

- y Qv^ • W>* + V{x)\iP\ 2 + |^| 4 j dx + 



i/=0 



v-0-v^ + y(x)iv;i 2 + 



-V^-V^* + y(x)|^| 2 + |V'| 2 |^| 2 )df-L I j / |^| 2 dx 



^) >o, 



(40) 



where we have used (|33|). This contradicts the quadratic decrease of Tt(ty(is)) near = 0. Thus, 
no ip exists which satisfies ||"0|| 2 = C and which has a lower energy value than ip. This proves the 
theorem. I 



Appendix B: Collapse and blow-up 

A comprehensive overview of many aspects of collapse and blow-up for the attractive (a = —1) 
nonlinear Schrodinger equation (|lj) with V(x) = is given in [33]. However, some modification 
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of the theory is required for nonzero potential V(x) 7^ 0. In p9, [Ml such modifications were 



considered for harmonic trap potentials. However, part of the analysis given in [34| applies to 
general potentials. This is the starting point of our considerations. 

For the remainder of this appendix, we consider condensates with attractive atomic interactions 
for which a = — 1. Equation (|l|) has two conserved quantities, given by the L 2 -norm of the solution 
(this has an interpretation as the number of particles per period N) and the Hamiltonian (|30|), In 
quantum mechanics, the nonlinear term in (|l|) is absent, and i/j(x,t) has the interpretation of a 
probability density function. This allows the definition of the expected value of a physical quantity. 
For instance, 

X =< x >= — I x \if)\ 2 dx (41) 

J fz 

is the expected value of the position. Such quantities can still be constructed in the presence of the 
nonlinearity, although some care is required to interpret them. The above definition leads to 

^ = - (W(x)) . (42) 

This equation is Newton's law. In quantum mechanics, the expected values of physical quantities 
satisfy classical equations of motion. That this equation persists in the presence of a nonlinear 
term is surprising. 



To examine collapse, the quantity |33| 




W ® = J n [f^ Xi ) H df (43) 
is studied. Note that it is by definition a positive quantity. The solution of (Q) is said to collapse 



if W(t) becomes negative [33, 34]. Use of the conserved quantities and Newton's law ( f42| ) gives the 
variance identities 

= 2 f \^\ 2 dx-d [ \^\ 4 dx-2 [ \?p\ 2 (x- VV{x))dx (44a) 
dt Jn Jn Jn 

4j- o / i ,12/ 



= 4W(V>)-(d-2) / |^| 4 dx-2/ \^\ 2 (2V(x) + x-VV(x))dx (44b) 
Jn Jn 

These identities are generalizations of variance identities for the free nonlinear Schrodinger equation 
(V(x) = 0). With d > 2, 

d 2 W f 

< mty) - 2 \^\ 2 G{x)dx, (45) 

with G{x) = 2V(x) + x ■ VV(x). For specific choices of the potential (such as a harmonic trap, for 
which Q = M. d ) a more detailed analysis of the collapse is possible. However, for the potential (|3|), 
the analysis here is restricted to general statements. 
If G{x) > then 

d 2 W 

-^r ^ (46) 

=> W(t) < 2H{^)t 2 + (3t + 7. (47) 

If for the initial condition TL{ip) < 0, then for all t, Tt(ip) < 0, since the Hamiltonian is conserved. 
Then d 2 W/dt 2 < 0, and for some time t c , W(t) < if t > t c . A crude overestimate of this time t c 
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is found from (|47|). Presumably it is possible for solutions with positive energy to collapse as well, 
but this requires more subtle arguments |m| , j35[ l than those used here. 

Thus, for a given potential one checks whether G(x) > 0. If so, then every initial condition of 
(|l]) with negative energy collapses in finite time. Note that it is always possible to ensure G(x) > 
on f2: this can be achieved by adding a constant Vq to the potential V(x), since is a bounded 
domain. However, this increases the energy of the solution by an amount NVq. Note that the 



absolute level of the potential is physically irrelevant, which is clear from ( 44a ), where only the 
gradient of the potential appears. In order to maximize the set of initial conditions for which 
collapse can be concluded, the constant Vq is chosen such that G(x) > 0, where equality is achieved 
in f2. 

The exact solutions constructed in Section || are stationary, thus W(t) is constant for these 
solutions and d 2 W/dt 2 = 0. Thus, no conclusion for these solutions follows from our analysis. 
By increasing the amplitude of the initial conditions, we no longer find exact stationary solutions. 
However, it is possible to show that these initial conditions lead to finite-time collapse. This is the 
approach followed for many of the numerical runs in two and three dimensions. 

Finally, following it is easy to show that collapse of a solution implies blow-up of that 
solution, i.e., sup^ \ip\ — ► oo. The argument here is easier than the argument in |34j because Q, is 
bounded. 

Theorem 4 For condensates with attractive atomic interactions (a = —1), collapse implies blow- 
up. 

Proof Consider 

2 / \ 2 

|2 



/ \^\ A dx< (suplV^ / 
Jn \ n / Jn 



dx = N sup |^| ) . (48) 
n \ n " 



The Heisenberg uncertainty principle states 



W(t) J VV" Vij*dx > ( — J , (49) 



thus 



2 , r 



sup |^| J > — dx 

£1 J J Q 



-j-(-2W(^)+ f W-V^*dx + 2 / V(x)\i;\ 2 dx 
N \ Jn Jn 



> -[-2H^)+l^pj 2 ± + 2 lj(xM 2 dx |. (50) 



Since J n V{x)\il)\ 2 dx G [— A^sup^ V(x), iVsup^ V(x)], this implies that as W(t) — > 0, sup^ l^l — > oo. 
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